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Abstract 

We have performed 2.5D and 3D simulations of conical jets driven by the rotation of an ordered, large-scale magnetic field in a 
stratified atmosphere. The simulations cover about three orders of magnitude in distance to capture the centrifugal acceleration as 
well as the evolution past the Alfven surface. We find that the jets develop kink instabilities, the characteristics of which depend on 
the velocity profile imposed at the base of the flow. The instabilities are especially pronounced with a rigid rotation profile, which 
induces a shearless magnetic field. The jet's expansion appears to be limiting the growth of Alfven mode instabilities. 
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1. Introduction 

Strong magnetic fields on large scales may play an essential, 
active role in the formation and evolution of jet-like outflows. 
The general idea is that a poloidal magnetic field, embedded in 
a plasma and anchored e.g. in an accretion disk or a black hole, 
is forced into rotation at the anchor point, a toroidal field de- 
velops and the plasma is accelerated by what can be interpreted 
as a centiifugal force in a corotating frame ( [Blandford & Payne 
1982) [Spruit 1996; |Konigl & Pud ritz 2000). However, this mag- 
netocentrifugal acceleration is only effective up until the Alfven 
surface, defined as the surface where the flow velocity equals 
the Alfven velocity. Beyond this point, the magnetic field will 
be strongly wound-up. Such a field configuration is potentially 
unstable with respect to certain MHD instabilities. 

MHD jets are susceptible to a variety of instabilities. Kelvin- 
Helmholtz (KH) instabilities are fed by the relative kinetic en- 
ergy between the jet and the ambient medium. They can dis- 
tort the jet surface only (ordinary modes) or the whole beam 
(e.g. Birkinshaw 1991] ), provoking shocks, mixing with ambi- 
ent material and possibly a disruption of the jet (Bod o et al.| 
1995 1998[ l. The presence of strong magnetic fields, poloidal 
or toroidal, is expected to hamper the growth of KH instabilities 
( |Appl & Camenzind|1992||Keppens et al.|1999| l. 

The free energy associated with the toroidal magnetic field 
is responsible for another class of instabilities, which is tradi- 
tionally known as current-driven (CD) and of notorious impor- 
tance in controlled fusion devices (for an introduction, see, e.g. 
[Freidberg 1987; Bateman 1978). The relevance for magnetized 



instabilities might destroy the ordered, symmetric state of a jet, 
leading to its disruption or, through magnetic reconnection, the 
associated dissipation of magnetic fields and steepening of the 



astrophysical jets has been pointed out by ,Eichler,(1993) ); Spruit 
|et al.| ( |1997| l; |Begelman| ( |1998| ) and others. Among CD instabiU- 
ties, m - 1 kink instabilities are the most effective. An ideal kink 
mode is characterized by helical displacements of the cylindrical 
cross sections of a plasma column. It is expected to grow on a 
dynamical time scale with respect to an Alfven wave crossing 
the unstable column. The susceptibility is strongly dependent on 
the magnetic pitch, a measure for the degree of wind-up. Kink 



magnetic pressure gradient, to its acceleration (Drenkhahn 2002 



Giannios & Spruit 2006). 
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Different kinds of instability can mix and interact. For exam- 
ple, Baty & Keppens ( 2002 1 show how CD instabilities can stabi- 
lize KH vortices at the jet surface. For this work, we used condi- 
tions under which CD kink instabilities are expected to dominate 
(low plasma-yS, small magnetic pitch). 

For a self-consistent study of kink instabilities, numerical 
simulations need to be carried out in 3D. ' Bell & Lucek| ( |1996| l 
did so using a simple model in which a toroidal magnetic 
field configuration was allowed to expand into a uniform atmo- 
sphere. This generated a jet which was subject to kink instabil- 
ities. [Nakamura & Meier] (|2004) performed 3D simulations of 
MHD jets in variously stratified atmospheres, finding that they 
can develop kink-like distoitions in the trans-Alfvenic region. 
Laboratory experiments of MHD jets have been performed by 
Hsu & Bellan ( 2005*), confirming that the magnetic pitch plays a 
crucial role for the formation of kink instabilities. 



1.1. Effects of jet expansion 

Jets from protostars, and especially AGN and microquasars, ex- 
pand in width d by orders of magnitude after passing through 
their Alfven radius. In an expanding flow there is no clean sep- 
aration between time dependence due to instability and that due 
to the expansion itself, making the question of stability less well 
defined. Analytical studies thus tend to focus on instabilities in 
a cylindrical geometry, with constant diameter (such as in the 
"magnetic tower" picture of Lynden-Bell 2003). Expansion has 
strong consequences on the behavior of instabilities, however, 
compared with jets modeled as cylinders of constant width. 

First, there is the tendency for the toroidal (azimuthal, 
around the jet axis) component of the magnetic field to dominate 
in an expanding jet. From the induction equation, the poloidal 
and toroidal components of the field vary as Bp ~ l/d^ and 

~ l/(i respectively (for constant jet velocity). Expansion thus 
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causes a continual increase of the ratio B^/Bp. Even when dissi- 
pation were to decrease the toroidal field at some point, the ratio 
increases again on further expansion. Free energy available in 
the toroidal field thus remains the dominant form of magnetic 
energy, and one may expect the question of stability and dissipa- 
tion to remain relevant on all length scales. It also follows that 
the nonlinear development of instabilities in an expending jet is 
expected to be very different from the constant-diameter case. 

Secondly, expansion has a strong effect even on the condi- 
tions for occurrence of instability. It has a stabilizing effect, since 
magnetic instabilities become ineffective when their signal speed 
(the Alfven speed) drops below the lateral expansion speed of 
the jet. This is discussed further below. 

1.2. Rationale of the calculations 

The aim of the calculations reported here is to study how kink 
instability operates under these conditions of expansion of the 
jet over several orders of magnitude in width. 

The degree of instability to be expected in a jet driven by 
a rotating magnetic field is intimately tied to the way it is 
collimated. If, instead of being cylindrical, the jet has a non- 
vanishing opening angle )?, the expected incidence of instability 
depends on the details of the dependence of § on distance r. An 
opening angle increasing with distance reduces instability, while 
for an asymptotically vanishing opening angle instability must 
always set in at some distance, if it was not present akeady from 
the start (see discussion in Sect. |1.3| l. 

The setup in the simulations presented here produces jets in 
the intermediate case of an (approximately) constant opening an- 
gle: a "conical" outflow. It turns out that in this case the presence 
of instabilities and their amplitude depends on secondary condi- 
tions such as the rotation profile imposed at the base of the flow, 
hence it is a good test case for the incidence of instabilities. 

Since the observed jets travel over such large distances, even 
marginal forms of instability can become effective. An important 
goal of the present calculations is therefore to cover a large range 
in distance, about 3 orders of magnitude. This is achieved by the 
use of a grid adapted to the approximately conical shape of the 
jet. 

1.3. Expected instability growth in expanding jets 

In the following, we estimate how the sideways expansion af- 
fects the growth of instabilities in broadening jets. In spherical 
coordinates (r, i?, if), we assume that the jet radius (distance to 
the jet's central axis) is given by 



where va,/) = B^/ yjAnp is the azimuthal Alfven velocity. The 
expansion rate is estimated by 



with R' - r sin ■&' 



(1) 



where the prime stands for evaluation at a reference distance r', 
for which we take a distance somewhat beyond the Alfven ra- 
dius. The flow has then approximately reached its asymptotic 
speed Vr ~ const, and the magnetic field has become predomi- 
nantly azimuthal. In the absence of magnetic dissipation due to 
instabilities, the field strength then varies as B^ cc (mag- 
netic flux conservation) and the density as p oc R^^ (mass con- 
servation). Since the growth rate Fg is expected to scale with the 
Alfven crossing rate v^^/R, we introduce a dimensionless insta- 
bility rate x of order unity: 



d\nR _ 1 drdR _ av^ 
df RAtdr r 



We find 



Tp va sin d-' 



AT 



(3) 



(4) 



with V :- Vi-Iva<p ~ const according to the ballistic approxi- 
mations mentioned above. Consequently, the instability growth 
rate dominates at some distance r for a collimating jet (a < 
1). Decollimation (a > 1), on the other hand, thwarts the 
growth of instabilities. A conical jet [a - 1) constitutes a lim- 
iting case where all depends on the combination of parameters 
x/(v sin j?')j which is of order unity. A numerical simulation is 
necessary to find out whether the instability or expansion pre- 
vails. 

The paper is organized as follows. In Sect. [2] we introduce 
the magnetocentrifugal jet model and account for the assump- 
tions made in our simulations. A detailed description of the nu- 
merical setup, the coordinate system and the scale-free units em- 
ployed in the analysis is given in Sect. |3] In Sect. [4] we give the 
parameters of the simulated cases and in Sect. |5]we present the 
results. There, we start by making predictions on the character- 
istics of instabilities by examining the relevant properties of our 
simulated jets. We proceed with an analysis of the instabilities 
that actually appeared and complete with looking for effects on 
the jets' dynamics. We finish with a discussion and conclusions 
in Sect.m 

2. The model 

The model is construed to apply to jets produced by ordered 
magnetic fields anchored in an accretion disk. This has become 
the default interpretation for the jets observed in AGN, mi- 
croquasars and protostellar objects, though it must be kept in 
mind that observational evidence of the key ingredient in this 



model, the ordered field ( 


Bisnovatyi-Kogan & Ruzmaikin|1976 


IBlandford & Payne|1982 


1, is still somewhat indirect. 



More uncertain is the shape of this field. The strength of the 
field anchored in the disk is likely to scale in some way with the 
orbital kinetic energy (or gas pressure) in the disk, hence will 
decline with distance R from the rotation axis. In the absence 
of more detailed information, we consider a simple form for a 
field of this kind, one in which the vertical (normal to the disk) 
component at the surface B, varies as B,{R) oc [1 h- (R/zq)^]^'^. 
Neglecting gas pressure and fluid motions, the field above the 
disk would be a potential field, its shape defined uniquely by B,. 
For V = 3/2 it is the field of a monopole with the source at a 
depth zo below the center of the disk. 

The initial state of the model is a gas distribution in hydro- 
static equilibrium in a field of this monopolar shape. Rotation is 



'~r' 



'' R' \r') 



(2) 



applied at the lower boundary, in a region R < Rq (see Sect. 3.3 
for details). This generates an outflow with an approximately 
constant opening angle on the order Rq/zq (a "conical" outflow). 
The surrounding volume remains approximately in static equi- 
librium, and serves to coUimate the outflow to the desired open- 
ing angle. 

The magnetic field responds to the rotation by winding up. 
That is, a toroidal magnetic field B^ is produced and gives, to- 
gether with the poloidal field Br, rise to helical field lines. The 
magnetic pressure gradient (minus the tension force) associated 
with B^ gives rise to a poloidal acceleration of the plasma which 
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is of centrifugal nature in a corotating frame. Beyond the Alfven 
radius, the acceleration ceases to be effective, while the field be- 
comes predominantly azimuthal. The further development de- 
pends on how strongly the jet is affected by instabilities in this 
highly wound-up field. Possibly, they endanger the jet's integrity 
and/or facilitate magnetic reconnection events. Magnetic field 
dissipation can entail further acceleration of the jet (Dre nkhahn] 
|2002) . As discussed above (Sect. 1.3 1, the "conical outflow" pro- 
duced in our monopolar background field is of special interest as 
it marks the boundary between cases expected to be strongly re- 
spectively weakly unstable. 

A self-consistent investigation of the problem requires a full 
3D treatment, because kink instabilities are non-axisymmetric. 
In addition to every 3D simulation we also performed an ax- 
isymmetric (2.5D) simulation with the same boundary and ini- 
tial conditions. This way, we could detect whether the jet evolves 
differently due to the instabilities. 

The basic parameters of the model are the magnitude of 
the rotation velocity, its profile Q(R), the relative strength of 
the magnetic field as measured by plasma-yS of the initial state, 
and the jet's opening angle. The parameter values are chosen 
such that the Alfven radius of the resulting outflow is located 
within the computational volume, so that the centrifugal acceler- 
ation process is covered in the simulation, but close to the lower 
boundary, so that the subsequent evolution can be followed over 
as large a distance as numerically feasible. Increasing the im- 
posed rotation rate moves the Alfven radius toward the lower 
boundary. Due to numerical limitations, however, could not 
be increased indefinitely in the simulations, and a compromise 
was necessary. In the results reported below the region inside 
the Alfven radius occupies about 10-20% of the box length. 



3. Methods 

We employ a spherical grid for our jet simulations. This enables 
us to follow jets with opening angles over a much longer dis- 
tance than is possible with a Cartesian grid, because the jet need 
not be overresolved at large heights in order to properly resolve 
its base. The obvious choice of letting the jet propagate along the 
polar axis is numerically problematic if non-axisymmetric flows 
are involved, because the grid is singular there. We therefore let 
the jet flow in equatorial direction. The computational volume 
covers a range AO = A<p in the polar and azimuthal angles, ad- 
justed to the opening angle of the flow. 

3.1. MHD equations and numerical solver 

We numerically solved the ideal adiabatic MHD equations, in- 
cluding a temperature-dependent temperature-control term K - 
K(T(t)). Explicitly, the equations are: 



5e 
5f 



5v 



+ V ■ 



-I- V ■ Vv 



dp 

^ + V ■ (pv) = 0, 
at 

--Vp + -^(V X B) X B - Vd) 
p A-np 



e + p + — \v B{B v) 



dB 

~dt 



= V X (v X B), 



where 



1 2^5' P 

e - -pv H h 

2 



Stt y-\ 



(5) 
(6) 

pvV<b + K, (7) 
(8) 

(9) 



is the total energy density, y - 5/3 is the adiabatic index, p is 
the gas pressure, O is the gravitational potential (external, no 
self-gravity) and the other symbols have their usual meanings. 
A notorious problem with low-/? MHD simulations in fully con- 
servative form, as in the code used here, is the amplification of 
discretization errors that occurs because the gas pressure is only 
a small contribution to the total energy (cf. |9]l, which is domi- 
nated by the magnetic energy. As in the case of highly supersonic 
flows, these errors manifest themselves in the form of "negative 
pressures" at occasional grid points. This problem does not oc- 
cur when an equation for the thermal energy equation is used 
instead of the total energy. We compute, in parallel, an alterna- 
tive update of the gas pressure from the thermal energy equation, 
in the form 



dp 

dt 



+ V ■ ipv) - -(y - l ypV v + K. 



(10) 



Where negative pressures appear, they are replaced by this value. 

Another device that turns out very useful to avoid negative 
pressures is the temperature-control term K in Eq. (|7]i. For this 
we use a scheme loosely modeled after Newtonian cooling, or an 
optically thin radiative loss process. After every full time step, 
we add/subtract thermal energy according to 



Ap{t) _ _ T{t) - Tit = 0) At 
Pit = 0) " Tit = 0) V^' 



(11) 



where T is the temperature and tk is a time scale chosen so as 
to keep the temperature within about a factor 30 of the initial 
atmospheric value. 

With the MHD Poynting vector 



Eq. (|7]i can also be written as 



ivxB)x B, 



5(e-Hp(D) 
di 



1 2^ r 

-pv + - 

2 T - 1 



p + p<^\v +S 



(12) 



= K, (13) 



describing the change of total energy including gravitational po- 
tential energy. We will employ this form later when we look at 
the energy flow rates. 

We used a newly developed Eulerian MHD code 
(Obergaulinger et al., in preparation) to solve Eqs. (|5]-[8]l. 
It is based on a flux-conservative finite-volume formulation 
of the MHD equations and the constraint transport scheme to 
maintain a divergence-free magnetic field ( [Evans & Hawley 
1988[). U sing high-resolution shock capturing methods (e.g. 



LeVeque 1992|), it employs various optional high-order re- 
construction algorithms and approximate Riemann solvers 



based on the multi-stage method ( |Toro & Titarev 20061. The 
simulations presented here were performed with a fifth order 
monotonicity-preserving reconstruction scheme (Suresh & 
Huynh 1997) , together with the HLL Riemann solver (jHarten 
[1983 1 and third order Runge-Kutta time stepping. 



3.2. Grid coordinates 

In the 3D simulations, our computational domain was centered 
around the y-axis in a "standard" spherical coordinate system 
(r, 6, 0) with 6 being the polar angle from the z-axis and <p being 
the azimuthal angle about the z-axis in the x-y-plane. See Fig.[T] 
for an illustration. Positioning the domain in equatorial (rather 
than polar) direction yields a grid which is free of singularities 
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Figure 1. Computational domain and coordinate nomenclature 
(schematic drawing). The jet propagates along the y-axis (9 - 
(p - njl), near which the grid is quasi-Cartesian in the normal 
plane. To describe the results of the simulations, we use the al- 
ternate spherical coordinate system (r, i?, ip) shown in the upper 
right inset, which takes the y-axis as the polar axis. 



and quasi-regular in transverse jet direction: Near the y-axis, the 
distance between two - const curves is Ax « r if we ne- 
glect terms of order {6 - njl)^ and of order (jp - nITf. The dis- 
tance between two 9 - const curves is then Az ~ r A9. With uni- 
form spacings A9 and A0, we thus obtain a grid whose area el- 
ements are approximately those of an equidistant Cartesian grid 
in a plane normal to the y-axis. 

For data analysis and plotting we use the "alternative" spher- 
ical coordinate system (r, i?, (p) with & being the polar angle from 
the y-axis and ip being the azimuthal angle about the y-axis in the 
z-x-plane. It is adapted to the propagation direction of the jet and 
as such better suited to describe its physics. R - rsm & denotes 
the distance to the y-axis. Generally, we refer to the r-direction 
with "poloidal" or "radial", the i/3-direction with "toroidal" or 
"azimuthal" and the y-axis as "polar axis" or "central axis". 

In the 2.5D simulations, the jet propagates along the z-axis 
and the ^-direction is taken to be symmetric. However, to avoid 
confusion, we use the same nomenclature as in the 3D simu- 
lations throughout this paper That is, we visualize the jet as 
propagating in y-direction and employ the (r, §, (p) system for 
describing its physics. 

For a proper resolution at all radii, it turned out to be neces- 
sary to employ logarithmic grid spacing (fPark & Petrosian 1996| l 
in r-direction. The r-left interface of grid cell i is situated at 



r ^n-l \'/" 



(14) 



where n is the total number of cells in the domain which is 



corresponding to a magnetic monopole of charge g located at the 
coordinate origin. The associated vector potential 



(16) 



was employed in the numerical setup to ensure solenoidality of 
the discretized magnetic field. Satisfying V x B = 0, the initial 
magnetic field is force-free. 

We impose the static gravitational field of a point mass M, 
also located at the origin. The stratification of gas pressure in 
this potential is chosen such that the plasma-/? := p/pmag in 
the initital state is constant throughout the computational do- 
main. Hence, since B oc r p oc r^*. The density in the ini- 
tial state is determined from hydrostatic equilibrium: pGM = 
-r^dp/dr oc r^^ where G is the gravitational constant. The tem- 
perature, sound speed and Alfven velocity vary as T oc r"', 
Cs oc r"'/2 and va oc r"'^^ in this stratification. 

The lower boundary of the computational volume, where the 
jet's "nozzle" resides, is at a distance r„ from the origin. The 
conditions at this surface are related to the gravitational potential 
by 

0(.) = -^ with M=^, (17) 
r Gpn 

where the subscript n (for "nozzle") denotes the values at rn. 

At the sides (9 and (p) and top (upper r) of the domain, we 
use open boundaries which allow for an almost force-free out- 
flow of material: p, p, all components of v and the transverse 
components of B are miiTored across the boundary interface to 
the opposing "ghost cells", the normal component of B is deter- 
mined by the solenoidality condition. The open boundaries work 
well and cause only minimal artefacts in the form of reflections. 
At the bottom (lower r) of the domain, where the jet emanates, p 
and p are kept fixed at their initial values in all ghost cells. The 
magnetic field is determined by extrapolation from the interior 
of the domain using the same scheme as for open boundaries. 
The velocity is prescribed to be zero except for the ghost cells 
below the nozzle area (R < Rjet^„ at r - r^). There, an azimuthal 
velocity field v = v^e^ is maintained, with either a Keplerian 
velocity profile 



j-^max for 0.2/;jet,„ <R< -Rjet.n 



elsewhere 



(18) 



or a rigid rotation profile 



lo 



jet,n 



elsewhere. 



(19) 



Note that we use the term "Keplerian" to indicate only that 
oc R^'^l^. The central mass M only serves to balance our cho- 
sen stratification and is not to be understood as the center of an 
accretion disk. 



bounded by r = and r - r"' 
cell / is situated at = (rj + r[)l2. 

3.3. Initial and boundary conditions 
The initial magnetic field is 



r" The cell center of grid 3.4. Units 



The setup described above is unambiguously determined by 6 



parameters, B^, p„, p„, R 



jet.m 



and 



but they are not all 



B(r) 



(15) 



independent. As units of length, pressure and density we use 
to = 27?jet,n, P() = Pn and po = Pa. The physical quantities ex- 
pressed in these units are listed in Table [T] Since these units 
are arbitrary, the number of independent parameters defining the 
problem reduces to 6 - 3 = 3. These are a plasma-/? value (which 
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Table 1. Normalization units 



Quantity 


Symbol(s) 


TTnit 


length 




/a 


gas pressure 


n 
F 




density 





0(\ 


velocity 


V 


Cso = ^jypolpa 


time 


t,T 


to - 'o/CsO 


energy density 


e 


Po 


energy flow rate 


6 


Poll/to 


force density 


F 


Polk) 


magnetic flux density 


B 


Bo = ^^npo 


current density 


j 


jo = Boc/(4nlo) 



determines an opening angle &o :- arcsin(/o/2rn), and a 
Mach number for the rotation: either v^^/cs^n or v™''/vA,n- 

The sound speed, Alfven velocity and escape velocity at r = 
r„, expressed in the normalization units, are Cs,n - Cso, 



VA,n 



and 




(20) 



(21) 



For the sake of clarity, we usually omit the normalization 
unit. For example, v = 5 would denote a velocity of 5cso, which 
is sonic Mach 5 at the jet nozzle. 

4. Cases studied 

In the following we present the results of two 3D simulations 
for two different rotation profiles imposed at the nozzle, the 
Keplerian and rigid rotation profiles given by (18 case K3) and 



( [19] case R3). These are compared with two 2.5D simulations 
with the same initial and boundary conditions (cases K2 and R2). 

In all cases, the initial state has a constant plasma-yS of 1 /9 
(B„ = 3), the maximum rotation velocity at the boundary is 
v™" - 0.33cs n = 0.1vA,-,n and the initial (half) opening angle 
is i?o = 5.7° (r„ - 5). This choice of parameters yields a jet 
with a magnetic pitch low enough to be unstable to kinks. At 
the same time, it avoids numerical problems found to arise with 
higher (supersonic) rotation velocities as a boundary condition. 
The "grid noise" in the 3D simulations (the grid is not axisym- 
metric in the rotation direction (p) turned out to be sufficient to 
excite instabilities, so we did not need to apply a perturbation by 
hand. For the temperature-control term, we used = 2. 

In the 3D simulations, we used 384 logarithmically spaced 
grid cells in radial direction and 96 uniformly spaced grid cells 
in each of the two angular directions. The physical extent of 
the simulated domain was 500 in the radial direction and 33.8° 
in each angular direction. The ratio between the maximum and 
minimum r is 505/5 = 101. The 2.5D simulations were per- 
formed with the same resolution in the radial direction and 64 
grid cells in the evolved angular direction which had an extent 
of 16.9°. The 3D simulations each ran for about one week (wall 
clock time) on 64 processors with MPI parallelization. 



5. Results 

The jets are initially accelerated mainly by gas pressure. This 
holds up to the sonic surface, which lies about halfway to the 
Alfven surface. Then, the Lorentz force becomes the dominant 
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thick lines: enclosing 95% of (?(^ 
thin lines: enclosing 75% of (§,„. 



100 



200 



30C 
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Figure 2. Jet border, defined as the polar angle within which a 
given percentage of energy flows (see Eq. 26 for the definition 
of £tot), as a function of distance. In all cases, most of the en- 
ergy flow is contained within the § ^ -&„ = 5.7° surface, where 
i?o is the initial opening angle. A minor but increasing amount 
of energy flows outside this angle. The energy density in the jet 
(especially the toroidal field) causes it to decollimate somewhat 
compared with the conical configuration in which it is embed- 
ded. 



-1 - " 




Figure 3. Selected magnetic field lines in the 3D simulation with 
a Keplerian velocity profile. The right-hand plot shows the entire 
domain (r - 5 . . . 505), the left-hand plot only the lower part 
up to r » 200. The color coding represents the strength of the 
azimuthal magnetic field. The magnetic field lines, which were 
purely radial to begin with, wind many times around the central 
axis, rendering it susceptible to kink instabilities. 



driving force. The Alfven radii are at r « 30 . . . 140, depend- 
ing on the simulation and the direction §: near the axis (small 
§), the poloidal field B,- is amplified and the Alfven radii are at 
larger distances than in the outer regions (large i?), where Br is 
attenuated. The opening angle of the jets increase somewhat with 
distance, but to a first approximation the flow can still be treated 
as conical, see Fig. |2] The jet front crosses the upper boundary 
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Figure 4. Density and temperature in a meridional slice in sim- 
ulation K3 (Keplerian rotation profile, 3D). The density is en- 
coded in intensity, with bright colors representing regions that 
are overdense with respect to the environment. The (square root 
of the) temperature is encoded in the hue, with blue meaning 
cold and red meaning hot with respect to the environment. The 
hoop stress associated with squeezes the plasma towards the 
central axis and creates an underdense cavity around the central 
part of the jet. At the boundary of this cavity the environment ex- 
erts the stress that confines the jet. The observed opening angle 
for a jet like this would be smaller than the width of the cavity. 



(r - 505) at t ^ 260 in the 2.5D simulations and at t ^ 330 
in the 3D simulations. The latter are subject to more numerical 
dissipation of kinetic energy, because the grid there is not sym- 
metric in azimuthal direction. This reduces the injected Poynting 
flux, see Fig. 1 1 Apart from that, 2.5D and 3D simulations give, 
for our purposes, comparable results. Fig.[3]shows how the mag- 
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Figures. Magnetic pitch as a function of distance in the 2.5D 
simulations. The magnetic pitch is also the smallest possible 
wavelength of an instability. The dependence of h on the polar 
angle § is stronger in the Keplerian case (left plot), suggesting 
higher stability. 
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Figure 6. The expected onset of instability depends on the ra- 
tio of the azimuthal Alfven speed and the expansion velocity. 
This ratio is shown here for the 2.5D simulations. The horizontal 
dotted lines are for two estimates of the condition under which 
growth is possible (see text). Below the respective line, expan- 
sion prevails and an instability cannot grow. 



netic field is wound up inside the jet in one of the simulations. A 
typical density and temperature distribution is shown in Fig.|4] 

5. 1 . Expected instabilities 

The observed instabilities can be compared with expectations 
from linear stability theory. To do this, we extract from the 
axisymmetric simulations the quantities that enter the stability 
conditions, and then compare the result with the evidence of 
non-axisymmetric instability in the corresponding 3D simula- 
tion. The available stability conditions apply only to steady or 
static configurations and have been derived either in the context 
of controlled fusion or cylindrical jets ( Appl et aL]|2000 Lery 
et al.|2000 ), hence the comparison can only be indicative. 

According to the Kruskal-Shafranov criterion, the longitudi- 
nal wavelength of an instability must be at least as high as the 
magnetic pitch on the unstable surface, defined to be the distance 
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: = 904.1 R2 : = 90 




Figure?. Alfven crossing times with t = tt in the 2.5D sim- 
ulations. In both cases, and with it the expected instability 
growth time increases with radius and distance from the central 
axis. The curves for r x 100 represent an underestimate, because 
the jet still accelerates at this distance. 



covered during one revolution of a helical field line about the 
central axis. Besides being the result of linear stability analyses 
in the context of controlled fusion, it can be derived heuristically 
from geometric arguments (Joh nson et al.||T958l l. Therefore, it 
should give a convenient scale also in cases for which it was 
not originally intended, like the expanding jets studied here. 
Deviations can be expected e.g. from the effect of one-ended 
line-tying, which in some cases has been found to lead to in- 
creased instability as opposed to a configuration without a free 



end ( |Furno et al.|2006l|Lapenta et al.|2006||Sun et al.|2008] l. 

For a conical jet, the magnetic pitch is 



h = IjtR 



(22) 



on the § = const surface. See Appendix [A] for a derivation of 
this expression. In the simulations, h decreases with r and settles 
to a constant value above the Alfven radius, see Fig.|5](compare 
also Fig. |3]l. The variation of the pitch with § depends on the 
kind of rotation imposed at the lower boundary. In the Keplerian 
case, the dependence is strong, with the asymptotic pitch being 
approximately 10 near the axis and 40 at the limb of the jet. In 
the rigid rotation case, h x 25 in all directions within the jet. 

The Alfven crossing time in a conical, unaccelerated jet, de- 
fined as the time it takes an azimuthal Alfven wave to orbit the 
central axis, is given by 



iR 



(23) 



where l - 2n for a full revolution, va^ - const is the azimuthal 
Alfven speed and vr - Vr sin § is the expansion velocity. In the 
simulations, va^ ~ const above the Alfven radius. This is as 
expected theoretically from conservation of mass and magnetic 
flux in a conically expanding, steady axisymmetric jet. is finite 
and physically meaningful only if the condition 



VAf_ 



> L 



(24) 



is satisfied. If it is not, the expansion takes place too fast for 
an Alfven wave to cross the jet and an Alfven mode instability 
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Figure 8. Radial component of the current density (V x B) in the 
3D simulations just before the jets reach the upper boundary. The 
rod in the middle of the jets is their central axis (y-axis). Helical 
distortions, characteristic al for kink instabilities, can be seen in 
the backward current (blue) in both cases. The amplitudes and 
wavelengths are significantly larger in the simulation with a rigid 
rotation profile (right-hand image). 



cannot grow. While the critical value of t is arguable, we note 
that causal contact across the jet by Alfven waves is only pos- 
sible if i > TT. The situation in our simulations is illustrated in 
Figs. |6] and |7] Instabilities can grow only slowly on magnetic 
surfaces with large §. Depending on the l needed for efficient 
growth, they may even be stalled due to the jet's expansion. In 
any case, instabilities grow most rapidly if they start at small r. In 
regions where the jet is accelerating (dvv/dr > 0) or decollimat- 
ing (d)?/dr > along a field line), the effective Alfven crossing 
time is underestimated by Eq. ( |23] l. The jet is then stabler than 
condition (|24]i suggests. 



5.2. Instabilities found in the simulation 

We observe non-axisymmetric, kink-like distortions in the mag- 
netic field and other quantities in both 3D simulations. They 
emerge near the Alfven radius, propagate with the flow and grow 
in amplitude along the way. It is convenient to look at the current 
density j = ^ V x B for a quantitative analysis. The radial com- 
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Figure 9. Position of the barycenter of the backward current 
(blue material in Fig. [8]l in simulation R3. The left-hand map 
shows the amplitude of the instabilities and the right-hand one 
its phase. The blank region is where the jet has not been yet, its 
border marks the jet front. Observers moving with the flow fol- 
low a time-position curve which is (approximately) parallel to 
the front, i.e. the instabilities are at rest with respect to a comov- 
ing frame. 




Figure 10. Amplitude of the perturbations (crosses, left-hand 
axis) as seen by an observer (dotted line, right-hand axis) located 
just behind the jet front (upper panel) and located 130 length 
units behind the jet front (lower panel) in simulation R3. The 
exponential fits (solid line, left-hand axis) yield the instability 
growth time. 



In the rigid rotation case CR3), the distortions attained large 
amplitudes of several degrees. Looking at in the r - const 
plane, we find that the whole jet is affected by the kink. The 
number of visible radial nodes is corresponding to wave- 
lengths on the order of 150, i.e. several times larger than the 
magnetic pitch. Owing to the distortions in the magnetic field, 
the axial current was perturbed as shown in Fig. [8] on the right- 
hand image. j7 helically twines around the central axis in rem- 
iniscence of "ideal" kink instabilities with an azimuthal mode 
number m- 1 . 

To analyze the unstable displacements, we determine the 
barycenter of the backward current j^. in the r - const plane, 
denoting its location with {d-j, ipj). The result, from which one 
can directly read off amplitudes and wavelengths, is shown in 
Fig. |9] The slope of the points of constant phase ipj in the r-t 
diagram corresponds with the flow velocity v,.. That is, the in- 
stabilities are at rest with respect to a comoving frame. We esti- 
mate the growth time in such a frame by introducing an observer 
moving with flow and measure d-j in doing so. We find strictly 
increasing, exponential growth if the observer is located just be- 
hind the jet front, see upper panel in Fig. 10 The exponential 
growth time t„ is generally on the order of the Alfven crossing 
times shown in Fig. [7] For observers which are farther behind 
the jet front, the amplitude does not follow a simple exponential 
increase, see lower panel in Fig. 10 for an example. Rather, it 
saturates and even declines in some cases. The reason for this is 
not clear. We cannot rule out the possibility that there is stabi- 
lizing feedback from the upper boundary. Considering that the 
flow is super- Alfvenic there, this seems unlikely though. 

In the Keplerian rotation case (K3), the jet also exhibits kink- 
like distortions, see left-hand image in Fig. |8] However, the per- 
turbation amplitudes are much smaller, with §j attaining peak 
values of about 1.4° directly behind the jet front and only about 
0.5° farther behind. Unlike in the rigid rotation case, only inner 
regions of the jet are affected by the kinks, the jet border is rela- 
tively unharmed. The wavelengths are on the order of 25-50, i.e. 
there are more radial nodes than in the rigid rotation case. Even 
for an observer traveling just behind the jet front, the amplitude 
is not strictly increasing, but saturates and tapers off after an ini- 
tial rise with a growth time of Tg ^ 50, i.e. also on the order of 
the relevant Alfven crossing times in Fig. |7] 



5.3. Impact on dynamics and energetics 

Instabilities release magnetic energy by transforming it into ki- 
netic energy, but for dissipation in the sense of magnetic re- 
connection sufficiently small length scales have to develop. In 
ideal MHD simulations like ours, such dissipation is present in 
the form of numerical diffusion due to the effect of interpola- 
tions across adjacent grid cells. The effect cannot be modeled by 
Ohmic resistivity but can be quantified through secondary effects 
like changes in the energy fluxes. 

We did not find conclusive evidence of magnetic field dissi- 
pation provoked by the instabilities in the 3D simulations. For 
example, the flow of nonradial magnetic field 



ponent j,- is related to and is as such characteristic for the dis- 
tortions in the magnetic field. In the unperturbed case, it is con- 
centrated about the central axis and along the outer boundary of 
the cavity illustrated in Fig.|4] with respectively opposite orien- 
tation. In our simulations, is directed in negative (/^-direction 
and the axial current, accordingly, in negative r-direction. We 
denote this backward current with j,T, so that j,- = + j^. 



J ^Bl + BlvrdA (25) 

r-const 

shows asymptotically the expected ~r behavior (as B^ ~ r"' and 
A ~ r^) with minor fluctuations but no decreasing trend when 
compared to the 2.5D simulations. 
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Figure 11. Energy flow rates through r - const in 3D (thick 
lines) and in the 2.5D (thin lines) simulations with a rigid rota- 
tion profile. There is no clear evidence of additional conversion 
of Poynting flux to kinetic energy. The energy flow through the 
lateral boundaries is virtually zero at all times. The situation is 
similar in the Keplerian case (K3 and K2). 
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Figure 12. Radial velocities in the 3D (thick lines) and in the 
2.5D (thin lines) simulations. The jets do not accelerate above 
the Alfven radius, despite instabilities. The value of v, near the 
central axis in simulation R3, though rising, does not give con- 
clusive evidence of acceleration, because the jet's axis is strongly 
distorted by the instabilities. 



It is helpful to look at the energy fluxes. From Eq. ( pjj ), the 
energy flow rate in poloidal direction is 



1 2 y - 1 



dA. (26) 



Decomposing the integral from left to right, we identify the flow 
rates of kinetic energy £kin, thermal enthalpy Sthsm, gravitational 
potential energy figrav and magnetic enthalpy fimag- We plotted 
these in Fig.[TT|for the simulations R3 and R2. The total energy 
flow £tot rises for small r due to the temperature-control term 
K in Eq. ([7]l. The conversion of Poynting flux to kinetic energy 
flux in the 3D simulation looks qualitatively the same as in the 
2.5D comparison simulation. In particular, there is no evidence 
of an additional conversion of £mag to £kin due to magnetic field 
dissipation. This agrees with the fact that there is no additional 



acceleration of the jets, see Fig. 12 



6. Discussion and conclusions 

We have simulated magnetocentrifugally driven, conical jets 
over a range in distance of 1000 times the initial jet radius, in 
both 3D and axisymmetric 2.5D. The calculations extend to a 
factor of about 5-10 beyond the Alfven surface. The 3D jets de- 
veloped non-axisymmetric instabilities of the kink kind. 

The violence of the instabilities depends on the rotation pro- 
file applied at the base. With a rigid rotation {ocR) profile, the per- 
turbations grow to much larger amplitudes than with a Keplerian 
(ocT?"''^) profile. We suspect that the reason for the differing be- 
havior lies in the magnetic shear, defined as the variation of 
the magnetic pitch with distance to the axis. In the rigid rota- 
tion case, there is virtually no shear as opposed to the Keplerian 
case, for which the pitch increases with distance from the axis, 
see Fig. [5] A shear-free configuration is expected to be unstable 
to non-resonant modes, whereas a configuration with increasing 
pitch is expected to be unstable to modes with a resonant sur- 
face inside the jet ( |Appl et al.|2000t|Lery et al.|2000| ). This fits 
well with what we observe in the simulations, viz. that the kink is 
confined inside the jet in the Keplerian case. Heuristically speak- 
ing, the differing behavior could be attributed to the fact that the 
outer (high ff) layers of the jet, which are more stable (higher 
magnetic pitch), damp internally arising modes in the Keplerian 
case. 

In both cases, the longitudinal wavelength of the instabili- 
ties is ~5 times larger than the value of the magnetic pitch near 
the axis. The relation is qualitatively consistent with the find- 
ings of Appl et al. (2000jl for a cylindrical jet. The growth time 
of the instabilities is on the order of the Alfven crossing time. 
The exact relation is difficult to determine, because the cross- 
ing time as well as the location of the resonant surface can only 
be estimated. As the azimuthal magnetic field strength and with 
it the azimuthal Alfven speed decrease past the Alfven surface, 
opposing parts of the jets become causally disconnected from 
each other. Thus, the jet expands too fast for Alfven mode in- 
stabilities to grow. The effect is amplified if the jet is diverging. 
Recollimation, on the other hand, should boost the growth of in- 
stabilities. 

As found in other studies, the conversion of magnetic en- 
thalpy (Poynting flux) to kinetic energy is fairly efficient, on the 
order 70%. Dissipation of magnetic fields by internal instabili- 
ties is expected to contribute additional acceleration of the flow 
( Drenkhahn|2002 1. The calculations do not show a clear signa- 
ture of this process. It seems that either the observed region is 
too small, and/or the numerical dissipation of magnetic fields 
is too weak. Also, from a macroscopic point of view, the in- 
stabilities were not violent enough to bring together fields with 
an antiparallel component, as is necessary for magnetic recon- 
nection to occur. Moreover, most of the magnetic enthalpy was 
already converted in the magnetocentrifugal acceleration pro- 
cess. Therefore, even if there was magnetic dissipation, the effect 
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would not be dramatic. Nevertheless, we found that the mag- 
netic field gets significantly distorted by the instabilities. This 
should facilitate magnetic field dissipation further downstream 
but it may be necessary to extend the calculations to larger dis- 
tances to see the effect. 

It is tempting to compare the instability-related structures in 
the simulations with structures in observed jets. The 3D jet struc- 
ture in Fig. |4] for example, is reminiscent of the semi-regular 
patterns seen in Ha images taken of outflows from young stel- 
lar objects (YSO) like HH 34 ( Reipurth et al.|2002| ). There are, 
however, other possible interpretations of the observed structure. 
The wiggles in YSO jets could also be the result of a precessing 
or orbitally moving source ( Masciadri & Raga,2002^ and refer- 
ences therein). The symmetric nature of structures often seen in 
jet and counterjet (e.g. HH 212 Wiseman et al.|200T| ) suggests a 
modulation of the outflow speed or mass flux originating at the 
source of the outflow rather than an instability developing further 
away. The irregularities caused by the instabilities studied here 
are possibly more important for internal magnetic energy release 
inside the jet than for major observable structures like the knots 
and wiggles in YSO jets, though they are likely to contribute to 
these at some level as well. 
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Figure A.l. Magnetic field line with pitch /i on a conical surface. 



Appendix A: Magnetic pitch for a conical jet 



The radial progress of the field Une depicted in Fig. |A.l| is deter- 
mined by 



dr 



air) 



(A.l) 



where R - r sin is the distance to the polar axis and a{r) is the 
unsigned slope. Assuming that Br oc r^^ and oc r"' due to 
magnetic flux conservation, we can write 



a(r) — flo 



(A.2) 



By integrating the resulting expression we obtain the radial dis- 
tance covered after one revolution: 



ro+h 



2k 



J"dr = J aoRod(f 

ro 



h - InanR, 



(A.3) 



Alternatively, we could also define a local magnetic pitch h by 
taking a = ao = const for the slope. The result is 



h 



- exp (2too sin §) - I. 



(A.4) 



The difference between h and h turned out to be insignificant 
in our analysis. This is understandable since h ^ h for small 
flo sin 1?. 



